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HIGHLIGHTS 


•  Thermal  shock  resistance  was  analyzed  for  direct  flame  fuel  cells. 

•  Mass/charge/heat  transport  processes  and  thermal  mechanical  stress  were  coupled. 

•  Effects  of  cell  structures  (anode  or  electrolyte  supported  SOFC)  were  studied. 

•  Effects  of  flame  temperature  non-uniformity  and  heating  process  were  studied. 
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A  detailed  two-dimensional  model  of  direct  flame  fuel  cell  (DFFC)  was  developed  by  considering  the 
coupling  effects  of  heterogeneous  chemical  and  electrochemical  reactions,  electrode  microstructure, 
transport  processes  of  mass,  charge  and  energy,  as  well  as  the  thermal  mechanical  stress.  The  stress 
distribution  was  simulated  at  different  heat-up  rates  which  represent  the  typical  DFFC  and  the  common 
solid  oxide  fuel  cell  (SOFC)  operation.  Transient  temperature  field  and  associated  thermal  stress  distri¬ 
butions  are  determined  and  analyzed  for  two  different  cell  structures.  The  failure  probability  of  the  fuel 
cell  is  defined  and  estimated  by  employing  the  Weibull  statistic.  The  model  is  demonstrated  to  be  a 
useful  tool  for  understanding  the  mechanical  stress  distribution  within  a  DFFC  cell  and  for  the  cell 
structure  design  and  optimization.  The  results  reveal  that  the  failure  probability  of  an  SOFC  cell  plate 
working  in  flame  conditions  may  be  6  orders  higher  than  that  in  the  common  SOFC  operation  conditions. 
The  anode-supported  SOFC  shows  better  thermal  shock  resistance  compared  with  the  electrolyte- 
supported  SOFC.  The  uniformity  of  the  flame  temperature  is  vital  in  the  DFFC  system  since  the  non- 
uniform  distribution  of  the  flame  temperature  greatly  increases  the  failure  probability. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Direct  flame  fuel  cell  (DFFC)  is  a  kind  of  novel  fuel  cell,  which 
combines  the  flame  and  the  solid  oxide  fuel  cell  (SOFC)  in  a  “no¬ 
chamber”  setup  [1].  The  flame  consumes  the  oxygen  at  the  anode 
side  and  provides  the  reacting  heat  for  keeping  the  SOFC  operating 
temperature.  The  cathode  is  exposed  to  the  air  to  maintain  the  gas 
concentration  gradient  between  the  different  electrodes  [1-6]. 
Compared  to  the  conventional  dual-chamber  and  single-chamber 
SOFC,  DFFC  has  the  advantages  in  fuel-flexibility,  sealing,  simple 
setup  and  rapid  start-up. 
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Thermal  stresses  are  the  main  factors  for  the  failure  of  solid  oxide 
fuel  cells,  especially  for  DFFC  operation.  During  the  operation  of 
DFFC,  the  fuel  cell  plate  is  heated  up  much  more  rapidly  by  the  flame 
compared  to  the  temperature  rising  rate  in  common  SOFC  operating 
conditions  and  this  may  lead  to  significant  thermal  stress  to  the  cell. 
In  addition,  the  non-uniform  temperature  distribution  of  fuel  cell 
generated  of  common  flame  can  also  lead  to  the  variation  of  fuel  cell 
temperature,  which  will  lead  to  the  thermal  stress  of  SOFC  as  well. 
Thus,  it  is  significantly  important  to  analyze  the  thermal  shock 
resistance  and  failure  possibilities  of  SOFC  for  operating  in  typical 
flame  conditions,  which  will  be  crucial  for  novel  DFFC  prototype 
design  and  performance  optimization.  Currently,  although 
numerous  studies  have  been  carried  out  on  thermal  stresses  for 
SOFC  experimentally  and  theoretically,  more  specific  models  should 
be  studied  in  order  to  quantify  and  predict  the  magnitude  of  the 
thermal  stress  in  SOFC  operating  in  flame  conditions. 


378 


Y.  Wang  et  al.  /  Journal  of  Power  Sources  255  (2014)  377-386 


In  this  study,  a  mechanistic  model  for  DFFC  by  coupling  the  mass 
transport  processes,  the  energy  conservation,  the  electrochemistry 
and  the  mechanical  stress  is  developed.  The  thermal  shock  resis¬ 
tance  of  the  DFFC  configuration  is  studied.  Transient  temperature 
field  and  associated  thermal  stresses  are  determined  and  compared 
for  two  DFFC  structures  which  are  based  on  the  anode-supported 
SOFC  (ASSOFC)  and  the  electrolyte-supported  SOFC  (ESSOFC) 
separately.  The  influences  of  the  flame  non-uniformity  on  DFFC 
performance  and  the  stress  field  are  studied.  The  Weibull  statistic  is 
employed  to  estimate  the  failure  probability  of  the  DFFC  during  the 
start-up  and  operation  period. 

2.  Model  development 

2.1.  Model  geometry 

The  geometry  of  the  model  exhibits  a  two  dimensional  (2D) 
planar  type  DFFC  as  seen  in  Fig.  1.  In  the  DFFC  setup,  the  fuel  cell 
was  located  above  the  flame  with  the  anode  facing  the  flame  front 
and  the  cathode  exposed  to  ambient  air.  Two  different  cell  struc¬ 
tures  (ASSOFC  and  ESSOFC)  are  studied  and  the  dimensions  are 
given  in  Table  1.  The  thicknesses  of  the  cell  layers  (ta,  te,  tc  in  Fig.  1 ) 
are  chosen  as  the  configurations  which  are  commonly  used  in 
published  literature  [7-9]. 


Table  1 

Dimensions  of  the  fuel  cells. 


ASSOFC 

ESSOFC 

1 

50  mm 

50  mm 

fa 

1000  pm 

50  pm 

fe 

20  pm 

150  pm 

tc 

60  pm 

50  pm 

Ionic  charge  balance  at  the  anode: 
Ql,an5act,an9(Vion,an-^el,an)  ,  ^  (  ^eff  yti/  \-n 

^  ion, an  v  v  ion, any  — 


CTPB  tpb 

-  -  I  i  _H 2_  ,  v  CCO 

-  (  *0,an,H2  5ulk  +  *0,an,CO  5ulk 

LH2  lCO 


Sact,an  I  CXP 


^aneF(veli 


an  ^ion,an  ^ref,an 


RT 


-exp  - 


(1  n:)neF^Vei  an  ^ion,an  ^ref,an) 


RT 


(1) 


Electronic  charge  balance  at  the  anode: 


2.2.  Model  assumptions 

The  main  assumptions  are  listed  as  the  following: 

(1)  The  electrochemical  reactions  spatially  occurred  along  elec¬ 
trode  thickness  within  the  electrode.  The  reaction  active  sites 
are  assumed  to  be  uniformly  distributed  in  each  electrode 
layer.  The  two  conducting  phases  (electronic  and  ionic)  are 
considered  to  be  continuous  and  homogeneous  in  each  layer. 

(2)  The  charge  transfer  reaction  is  assumed  to  take  place  at  the 
three  phase  boundary  (TPB). 

(3)  A  mixture  of  FI2,  CO,  H20,  C02  gases  with  a  typical  flame 
composition  were  chosen  as  the  anode  gases  and  the  cathode 
gas  was  chosen  as  air.  All  of  the  gasses  are  considered  as  ideal 
gases. 

(4)  The  convection  flux  and  pressure  gradient  in  the  porous 
electrodes  are  ignored. 

(5)  Radiative  heat  transfer  is  neglected  within  the  whole  cell. 

(6)  The  temperature  load  is  considered  to  be  the  only  one  in 
stress  analysis. 

2.3.  Governing  equations 
2.3.1.  Charge  conservation 

The  governing  equations  for  SOFC  electrode  transient  charge 
balance  were  shown  as  follows  [10,11]: 


Ql,an*^act,and(^el,an  ^ion,an) 


at 


+  V-(  ^efan^^efan)  —  Qel,a 


—  Qio 


(2) 


Ionic  charge  balance  at  the  cathode: 

Ql,ca^act,caa(Vion  ca  -  Vel  ca)  _  ^  f  ^eff  \  _  n 

Qf-  a  ion,ca v  Mon,cay  —  Clio 


CTPB  / 

=  'O.caWa-^jk  exp 

c02  V 


fccneF(ve  u 


ca  ^ion,ca  ^ref,ca 


RT 


-exp  - 


(1  Ol)neF(Ve  jca  ^ion,ca  ^ref,ca) 


RT 


(3) 


Electronic  charge  balance  at  the  cathode: 

Ql,ca^act,caa(Vei  ca  —  Vi0n,ca) 


at 


V-(-<caVV/el,ca)  =  Qel.ca 


—  Qion,c 


(4) 


Air 
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Fuel-rich  flame 


mmmir 

Hydrocarbon/Air 
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Fig.  1.  DFFC  setup  (left)  and  simplified  2D  model  geometry  of  planar  SOFC  (right). 
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where  t  is  time,  Cdi  is  the  specific  interface  double-layer  capacitance 
between  electronic  and  ionic  conductor  phases;  Q.  is  the  transfer 
current  source;  Vei  and  Vion  denote  the  electric  potential  of  the  two 
conductor  phases  [10]. 

2.3.2.  Chemical  reaction  model 

Reversible  water-gas  shift  reaction  (WGSR,  Eq.  (5))  and 
reversible  direct  internal  reforming  reaction  (DIR,  Eq.  (6))  have 
been  taken  into  consideration,  which  is  called  internal  reforming  as 
well. 

C0  +  H20^C02  +  H2  (5) 


2.3.3.  Mass  balance 

It  is  generally  agreed  that  one  of  the  most  convenient  approach 
for  modeling  the  gas  transport  in  the  porous  electrode  is  the  dusty- 
gas  model  (DGM). 

The  DGM  could  be  formulated  as  follows  while  neglecting  the 
pressure  gradient  in  the  porous  electrode  and  the  implicit  rela¬ 
tionship  between  the  species  mole  fraction  and  molar  diffusion 
fluxes  was  provided, 


-cvx,-  =  £; 


' XjNi-XfNj \  Nt 


Df 


Dfff. 

kn,i 


i  =  1,2,3..., n 


(14) 


CH4  +  H2O^CO  +  3H2  (6) 

The  above  chemical  reactions  are  assumed  to  be  limited  in  the 
anode  which  can  provide  Ni  catalyst.  Haberman  and  Young’s 
model  is  adopted  to  describe  the  rate  of  reversible  WGSR  and  DIR 
[12].  The  expressions  of  reversible  WGSR  rate  (RShift>  mol  m”3  s_1) 
and  reversible  DIR  rate  (Rreform,  mol  m-3  s-1)  are  shown  as 
following. 


where,  the  effective  molecular  diffusion  coefficients  D?ff  and  the 
effective  Knudsen  diffusion  coefficients  can  be  calculated  by 
considering  the  porous  material  property  as, 


Df 


oooioit175(4  +  ^ 


\l/2 


V?/3  +  yl/3 


(15) 


Rshift  =  Kf  (ph20PC0  - 


(7) 


neff 

^kn,! 


(16) 


ks f  =  0.0171  exp ) >  m°l  m  3  Pa  2  s  1 

Kps  =  exp(  -  0.2935Z3  +  0.6351Z2  +  4.1788 Z  +  0.3169) 


Z  = 


1000 

T(K) 


^reform 


=  K 


F  (pch4Ph2o 


Ph2Pco\ 

Kpr  } 


(8) 

(9) 

(10) 

(11) 


where,  e  is  the  electrode  porosity,  t  is  the  tortuosity  factor,  V  is  the 
diffusion  volume  and  r  denotes  the  average  pore  radius. 

Eq.  14  can  be  formulated  as  the  explicit  relationship  between 
species  model  fractions  and  molar  fluxes  using  matrix  notation,  and 
eventually  we  can  get, 

^'^  +  vf-pw,^DyVxJj  =  R,  (17) 

The  equations  for  electrode  mass  balances  can  be  summarized 
as  follows  (Table  2): 


fcrf  =  2395  exp^  , mol  m3  Pa~2  (12) 

Kps  =  1.0267  x  1010  x  exp(  -  0.2513Z4  +  0.3665Z3 

+  0.5810Z2  -27.134Z  + 3.277)  (13) 

The  heat  generation/consumption  from  chemical  reaction  is 
considered.  The  values  are  given  in  Table  4. 


2.3.4.  Energy  conservation 

The  governing  equations  for  energy  balance  were  shown  as 
follows: 

^H  +  V(-Aefrvr)  =  Qheat  (18) 

where  Aeff  is  the  effective  thermal  conductivity  of  the  porous 
electrodes, 


Table  2 

Governing  equations  for  mass  conservation. 


Domain 

Anode 


Cathode 


Governing  equations 

8(W,WM) 4  v-(-panwH20  E  (DHi0jVXj))  =  MH2o -  3Rreform  - Rshfft) 
+  V  •  (-fcnwH2  E  (Dh2  jVXj))  =  IMh2  (3Rreform  +  Rshift  - 

a(8myco;)+v.(_panWc0z  (DCQ2 jVXj))  =  Mco;  preform  +  Rshift  + 

j= 1 

a(C“^WC0)  +  V '  (-ft,, „WC0  E  (Dcoj VXjO)  =  Mco  (Rreform  -  Rshift  - 

8(g“,’a1W°2)  +  V-(-RcaW0;EjUl(Do2jV^))  =  
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Aeff  =  — 2AS+  g  1  (19) 

2AS+Ag  +  3AS 

where,  \s  and  lg  are  thermal  conductivities  of  solid  phase  and  gas 
phase  materials  in  porous  electrodes.  Qheat  is  the  heat  source  term, 
and  can  be  calculated  as  the  following  equation  in  different  solid 
structures, 


I  Qohm- electrolyte 

\  Qohm  +  Ore v  +  Qi rr-  porous  electrode 


(20) 


where,  the  Qohm  is  the  ohmic  heat  and  can  be  calculated  as, 


Qoh 


m 


(^,  electrolyte 

j2  j2 

+  Oei’ Porous  electrode 


(21) 


where,  i  is  the  local  electronic  or  ionic  current  density,  a  is  the 
electric  conductivity. 

Qrev  is  the  reversible  heat  effects  of  entropy  change  in  electro¬ 
chemical  reactions,  and  can  be  calculated  as, 


Table  4 

Pore  structure  parameters  in  the  porous  electrode  [17]. 


Cell  layer 

Porosity 

Mean  pore 
radius  (pm) 

Stpb 

(m2  m"3) 

SlMi 

(m2  m-3) 

Anode 

0.36 

0.193 

2.14  x  105 

3.82  x  106 

Cathode 

0.36 

0.193 

2.14  x  105 

- 

variables  at  the  boundary,  respectively.  wg  an  is  the  mass  fraction  of 
the  flame  species  and  wgiCa  is  the  mass  fraction  of  gas  species  in  the 
cathode.  The  difference  between  Vceii,an  and  VCeii,ca  is  the  cell 
operation  voltage  in  the  calculation,  and  VCeii,ca  was  set  as  zero.  Tan 
in  the  table  represents  the  flame  temperature  and  the  convective 
heat  flux  boundary  means  n  (-AvT)  =  0.  The  mechanical  bound¬ 
ary  condition  “Free”  means  that  the  fuel  cell  is  not  constrained  and 
free  to  bend. 


2.5.  Model  parameters 

2.5 A.  Model  parameters  for  the  thermo- electrochemical  model 
In  Section  2.3.1  and  2.3.3,  the  effective  reaction  area,  TPB  area 
and  TPB  length  are  needed  in  the  calculation.  The  parameter  Stpb 


Qrev 


l'trans,ca,H2oSact,ca^  (Sh2  0 +2Se-  —  Sh2  —  SQ2-  )+l'trans, ca,C02Sact,ca7’ (^C02  +2Se-  — Sco- SQ2-  ) 

2F 


+  ^shiftQshift  T-  ^reform  Qreforrm  Cathode 


j'trans.anSact.arT (Sq2-  ^So2  2Se- ) 

2F  ’ 


anode 


(22) 


where,  S  is  the  molar  entropies  of  species.  Then  the  cathode  en¬ 
tropy  change  can  be  calculated  by  subtracting  anode  entropy 
change  from  whole  the  entropy  change  of  electrochemical  reaction. 
Qirr  is  the  irreversible  heat  generation  due  to  activation  polariza¬ 
tions  which  can  be  formulated  as, 

Qirr  “  I^Qell  =  I^Qionl  (23) 

2.3.5.  Mechanical  model 

Thermal  stress  in  the  three-layered  structure  (anode,  electrolyte 
and  cathode)  was  calculated  using  the  solid  model  within  Struc¬ 
tural  Mechanics  Module  in  the  finite  element  commercial  software 
COMSOL  MULTIPHSICS®.  It  should  be  noted  that  the  mechanical 
model  considers  only  the  thermal  load  and  the  reference  temper¬ 
ature  is  assumed  to  be  1283  K  following  Atkinson  and  Selcuk  [13]. 

2.4.  Boundary  conditions 

The  boundaries  of  the  model  geometry  are  labeled  in  Fig.  1.  The 
detailed  settings  of  boundary  conditions  for  the  partial  differential 
equations  are  listed  in  Table  3. 

The  boundary  conditions  “Insulation”  and  “Continuity”  mean 
that  the  partial  derivative  is  zero  and  the  flux  is  continuous  of  the 


Table  3 

Boundary  conditions. 


Boundary 

Ionic 

charge 

Electronic 

charge 

Mass 

balance 

Energy 

conservation 

Mechanical 

df^an/ac 

Insulation 

Ycell.an 

Wg,an 

Tan 

Free 

0Qelec/an 

Continuity 

Insulation 

Insulation 

Continuity 

Continuity 

Qf^elec/ca 

Continuity 

Insulation 

Insulation 

Continuity 

Continuity 

0QCa/cc 

Insulation 

Ycell.ca 

Wg,ca 

Convective 
heat  flux 

Free 

Others 

Insulation 

Insulation 

Insulation 

Insulation 

Free 

can  be  formulated  by  using  the  particle  coordination  number  in 
binary  random  packing  of  spheres  together  with  percolation  theory 
as  [14,15].  The  physical  properties  of  gas  mixture,  e.g.  specific  heat 
capacity,  heat  conductivity,  viscosity  et  al.  were  calculated  in  the 
well-correlated  equations  with  temperature  in  the  works  of  Young 
et  al.  [16].  The  typical  SOFC  cell  materials  (i.e.  Ni-YSZ/YSZ/LSM) 
have  been  considered  for  this  study.  To  be  more  specific,  the  ther¬ 
mophysical  properties  of  zirconia  stabilized  with  8%  mol  of  Y2O3 


Table  5 

Model  parameters  for  the  thermo-electrochemical  model  [11,12,18,19]. 


Parameters 

Value 

Electric  conductivity  (Sm  ’) 

Ionic  conductor  YSZ,  aYsz 

3.34e4exp(  - 1 0,300/T) 

Electronic  conductor  LSM,  <7Lsm 

4.2e7/Texp(-l  150/7) 

Electronic  conductor  Ni,  (rNi 

3.27e6-l  065.37 

Heat  conductivity  ( W  m  1  I<  1 ) 

Anode,  Aan 

6.23 

Cathode,  Aca 

9.6 

Electrolyte,  Aeiec 

2.7 

Density  (kg  m  3) 

Anode,  pan 

3030 

Cathode,  pca 

3310 

Electrolyte,  peiec 

5160 

Heat  capacity  (J  kg  1  K  1 ) 

Anode,  Cp,an 

595 

Cathode,  CPiCa 

573 

Electrolyte,  Cp>eiec 

606 

Parameters  of  chemical  reaction 

DIR  reaction  heat  Qreform,  J  mol  1 

-20,600 

WGSR  reaction  heat  Qshift,  J  mol-1 

41,000 

Parameters  of  electrochemical  reaction  kinetics 

Interface  conductivity  /?an//Jc a,  ^  1  m  2 

6.8el2/5.8el0 

Activation  energy  £act,an/£act,ca.  J  mol-1 

120,000/130,000 

Microstructures 

Anode  porosity/tortuosity/pore  radius,  pm 

0.36/5/0.193 

Cathode  porosity/tortuosity/pore  radius,  pm 

0.36/3/0.193 
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Table  6 

Elastic  properties  of  each  cell  material. 


Temp 

Anode  Ni/YSZ 

Electrolyte  YSZ 

Cathode  LSM 

E  (GPa) 

298  K 

76.75  [21] 

190  [22] 

41  [20] 

1073  I< 

63.4  [21] 

157  [22] 

43.4  [20] 

V 

298  K 

0.283  [21] 

0.308  [22] 

0.28  [20] 

1073  K 

0.286  [21] 

0.313  [22] 

0.28  [20] 

a  (xlO-5) 

298  K 

11.7  [24] 

7.6  [23] 

9.8  [23] 

1273  I< 

12.5  [24] 

10.5  [23] 

11.8  [23] 

(YSZ),  Lao.8Sr0.2Mn03  (LSM)  and  a  Ni-YSZ  cermet  which  is  reduced 
from  an  oxidized  form  NiO-YSZ  (50%  wt  of  NiO)  were  used  in  the 
model,  as  shown  in  Table  5.  These  parameter  values  were  based  on 
the  values  from  published  literature  and  the  temperature  depen¬ 
dence  of  the  heat  conductivity,  density  and  heat  capacity  of  all  solid 
materials  were  all  considered  negligible  in  this  model. 

2.5.2.  Thermo -elastic  material  properties 

The  elastic  modulus  and  Poisson’s  ratio  of  the  porous  cathode 
and  anode  have  been  shown  to  be  strongly  dependent  on  porosity 
[20,21].  It  should  be  noted  that  the  porosity  of  the  anode  heavily 
depends  on  the  oxidized  or  reduced  state  of  the  nickel  phase  and 
the  anode  is  considered  to  be  in  the  reduced  state  with  the  porosity 
of  36%  in  this  model.  Besides,  the  thermo-elastic  material  proper¬ 
ties  are  dependent  on  the  temperature  and  are  summarized  in 
Table  6.  A  linear  interpolation  is  assumed  between  the  data  points 
quoted  from  the  literature. 

2.6.  Model  solution 

The  calculations  were  performed  using  the  finite  element 
commercial  software  COMSOL  MULTIPHSICS®.  The  symmetry 
mapping  structured  quadrilateral  mesh  was  used  in  the  solution. 
The  mesh  consists  of  7150  and  9750  elements  separately  for  the 
ASSOFC  and  ESSOFC.  The  default  stationary  nonlinear  solver  (uses 
damped  Newton  method  with  direct  linear  system  solver  UMF- 
PACK)  of  COMSOL  MULTIPHSICS®  was  used  to  solve  the  thermo¬ 
electrochemical  model.  The  time  dependent  segregated  solver 
was  used  to  solve  the  time-dependent  thermo-mechanical  model. 


different  temperatures.  In  order  to  validate  the  mechanical  model, 
the  stress  fields  in  the  half-cell  consisting  of  a  porous  Ni/8YSZ 
anode  layer  and  a  dense  8YSZ  electrolyte  layer  were  firstly  calcu¬ 
lated  and  compared  with  the  experimental  data  in  Table  7.  It  can  be 
seen  that  between  room  temperature  and  1073  K  a  decrease 
from  -300  MPa  to  -50  MPa  of  residual  thermal  stresses  has  been 
measured  for  10-pm  YSZ  electrolyte  on  510-pm  Ni-YSZ  anode  [27], 
which  is  from  -384  MPa  to  -81  MPa  in  simulation.  Fischer  et  al. 
[28]  reported  a  compressive  thermal  stress  of  520  MPa  for  a  half¬ 
cell  with  an  anode  thickness  of  1.5  mm  and  an  electrolyte  thick¬ 
ness  of  10  pm,  which  is  higher  than  the  simulation  result  423  MPa. 
There  exists  a  difference  between  the  measured  thermal  stress  data 
in  experiment  and  our  simulation  results,  but  they  are  within  the 
same  order  of  magnitude.  The  difference  is  possibly  due  to  the 
difference  in  anode  composition,  anode  density,  and  electrolyte 
density  which  will  influence  the  Young’s  modulus  of  these  mate¬ 
rials  and  the  thermal  stresses.  It  should  be  noted  that  it  is  difficult  to 
directly  compare  the  experimental  data  and  the  simulation  results 
since  the  exact  geometrical  parameters  and  material  properties  are 
not  available  for  the  stresses  reported. 


3.1.2.  Validation  by  the  analytical  model 

To  further  validate  the  mechanical  model  in  a  whole  SOFC  cell 
plate,  an  analytical  model  developed  by  Hsueh  [26]  was  applied  in 
this  section.  In  the  analytical  model,  the  stress  distribution  in  the 
nth  layer  of  a  multi-layer  can  be  calculated  by  considering  the  other 
n  -  1  layers  to  be  a  single,  composite  bar  with  composite  material 
properties  derived  by  averaging  the  material  properties  of  each 
layer,  weighted  by  layer  thickness. 

A  multi-layer  strip  is  shown  schematically  in  Fig.  2,  in  which  n 
layers  of  individual  thicknesses,  q,  are  laminated  sequentially  on  a 
substrate  with  a  thickness,  ts.  The  subscript  i  denotes  the  layer 
number  from  1  to  n. 

The  relation  between  hi  and  q  is  described  by 

ht  =  Y)  tj  (f=l  ton)  (24) 

j=  1 

The  total  strain  in  the  system  is  formulated  as 


3.  Results  and  discussions 

3.1.  Mechanical  model  validation 

The  validity  of  the  proposed  thermo-electrochemical  model  in 
this  paper  has  been  validated  in  our  previous  studies  [25].  The 
mechanical  part  of  the  model  is  further  validated  in  this  paper.  The 
mechanical  model  is  checked  by  experimental  data  and  an 
analytical  model  developed  by  Hsueh  [26]  separately. 

3.1.1.  Validation  by  experimental  data 

In  the  experiments  of  Refs.  [27]  and  [28],  the  residual  thermal 
stresses  in  the  electrolyte  layer  of  a  half-cell  were  measured  at 


Table  7 

Comparison  of  simulation  and  experimental  residual  thermal  stresses  in  the  elec¬ 
trolyte  layer. 


Temp 

Electrolyte 
thickness  (pm) 

Anode 

thickness  (pm) 

Electrolyte  residual 
stress  (MPa) 

Experiment 

Simulation 

298  K 

10 

510 

-300  [27] 

-384 

Experiment 

Simulation 

1073  K 

10 

510 

-50  [27] 

-81 

Experiment 

Simulation 

298  K 

10 

1500 

-520  [28] 

-423 

e  =  c  +  Z  (for  -ts<z<hn)  (25) 

where  c  is  the  uniform  strain  component  given  by  (28),  tb  dictates 
the  location  of  the  bending  axis  (29),  r  is  the  radius  of  curvature  of 
the  system  (30).  It  should  be  noted  that  the  bending  axis  in  Hsueh 
model  is  defined  as  the  line  in  the  cross-section  of  the  system 
where  the  bending  strain  is  zero,  which  is  different  from  the  con¬ 
ventional  neutral  axis  which  is  defined  as  the  axis  in  the  cross- 
section  of  the  system  where  there  are  no  longitudinal  stresses  or 
strains. 


i 


Layer  n 

tn 

Layer  i 

ti 

Layer  1 

tl 

Substrate 

ts 

z=hn 

z=h„_i 


Z=h, 

z= 0 

Z=-ts 


Fig.  2.  Definition  of  the  variables  in  analytical  multi-layer  thermal  stress  model. 
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C  = 


fb  = 


Estsots  +  Y^=  l  Ejtfaj  )  A T 


Ests  +  1la  =  i  Ej tj 

-Estj  +  J2i  =  i  Ejtj(2hj_ i  +  tj) 
2(Ests  +  E/  =  i^ 


(26) 

(27) 


£sfs2(2ts+3tb)+E;Li£!ti 

[6  hj  i  +  6/i  j  i  tj + 2 1?  —  3 (2h  j| + tj )  ] 

3 

£s(c-asAT)ts2- 

EiLi£ifi(c-aIAT)(2hi_1+ti) 

(28) 


The  normal  stresses  in  the  layers  are  related  to  strains  by 
o i  =  E{ (e  —  ajAT)  (for  z  =  l  to  n,  i  =  s)  (29) 

If  the  system  has  a  planar  geometry  rather  than  a  strip,  as  is  the 
case  with  the  SOFC  model  the  biaxial  Young’s  modulus  E/  should  be 
used  instead  of  Ej 


where  v  is  the  Poisson’s  ratio. 

The  stresses  were  calculated  for  the  ASSOFC  and  DSSOFC  using 
both  the  analytical  model  and  the  mechanical  model  used  in  our 
study  when  it  is  free  to  bend  under  the  application  of  a  uniform 
temperature  increase  of  1000  K.  In  the  analytical  model,  the  anode 
is  taken  as  the  substrate  layer,  the  electrolyte  as  layer  1  and  the 
cathode  as  layer  2.  The  elastic  properties  are  taken  to  be  uniform 
and  constant  throughout  a  layer  according  to  Table  6.  The  com¬ 
parison  of  the  two  different  solutions  is  shown  in  Fig.  3.  The  stress 
calculated  by  the  model  we  used  agrees  well  with  the  analytical 
result  with  a  difference  of  less  than  4%.  The  difference  is  due  to  the 
edge  effects  in  the  numerical  model  providing  stress  relief  but  the 
analytical  model  assumes  the  domain  to  be  infinite  in  the  in-plane 
dimension  [29]. 


3.2.  Stress  comparison  for  DFFC  configuration  and  common  SOFC 

The  DFFC  can  be  heated  up  to  the  working  temperature  within 
seconds  [5],  which  is  extremely  fast  compared  to  the  heating  time 
in  common  SOFC  operating  conditions.  In  this  section,  the  stress 
distribution  was  simulated  at  different  heat-up  rates  which 


Fig.  4.  Temperature  distribution  of  three  different  heat-up  rates. 


y  (f*m) 

Fig.  5.  Stress  field  distribution  of  DFFC  with  different  heat-up  rates. 


represents  the  typical  DFFC  and  the  common  SOFC  operation.  The 
heat-up  time  here  means  the  time  it  takes  the  anode  surface 
temperature  to  increase  from  room  temperature  (298  K)  to  oper¬ 
ation  temperature  (1073 1<).  Since  the  SOFC  is  heated  by  the  flame  in 
DFFC  condition,  the  heat-up  time  is  the  result  of  the  rate  of  the 
homogeneous  oxidation  reaction  and  the  thermal  mass  of  the 
system.  The  homogeneous  oxidation  reaction  is  rather  rapid  with  a 
characteristic  time  scale  of  several  ms,  and  we  considered  this  time 
scale  should  be  the  boundary  limits  (worst  situation  for  DFFC 
operation)  for  the  heat-up  time  of  the  anode  gases.  Flowever,  in 
practice,  the  heat-up  time  can  also  reach  the  time  scale  of  several 
seconds  ascribing  to  the  heat  transfer  resistance,  for  instance,  the 
convective  heat  transfer  between  the  flame  gas  with  the  cell 


Fig.  3.  Comparison  of  analytical  solution  with  mechanical  model  analysis  (a)  ASSOFC  (b)  ESSOFC. 
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Table  8 

Fracture  strength  properties  of  each  cell  material  [33]. 


Anode 

Electrolyte 

Cathode 

(Ni/YSZ) 

(YSZ) 

(LSM) 

Characteristic 

79 

154 

75 

Strength  (MPa) 
Weibull  modulus 

7 

8.6 

3.7 

surface.  As  a  result,  in  our  paper,  four  different  heat-up  rates  (0.1  s, 
1  s,  50  s  and  1000  s  from  298  K  to  1073  I<)  were  set  as  the  tem¬ 
perature  boundary  conditions.  The  former  three  conditions  repre¬ 
sent  the  time  scales  for  typical  DFFC  heat-up  processes,  which  is 
referred  as  “DFFC  case”  in  following  text  and  the  latter  one  repre¬ 
sents  the  time  scale  for  common  SOFC  heat-up  processes  is  referred 
as  “SOFC  case”  in  following  text.  It  should  be  noted  that  the  anode- 
supported  SOFC  configuration  is  adopted  in  this  section.  The  tran¬ 
sient  temperature  distribution  response  and  corresponding  stress 
field  when  the  anode  temperature  reaches  1073  I<  are  shown  in 
Fig.  4  and  Fig.  5,  respectively.  In  the  DFFC  case,  when  the  anode 
temperature  reaches  1073  K,  the  maximum  transient  temperature 
difference  between  the  anode  surface  (y  =  0)  and  the  cathode 
surface  (y  =  1150)  can  even  reach  several  hundred  degrees  (0.1  s 
and  1  s  case);  while  the  maximum  transient  temperature  difference 
across  the  cell  is  only  kept  within  1°  for  the  SOFC  operating  case. 
Flowever,  it  should  be  noted  that  the  temperature  difference  be¬ 
tween  the  electrode  surfaces  in  DFFC  case  can  be  controlled  within 
several  degrees  when  the  heat-up  rate  is  controlled  to  be  relatively 
slow  (about  1  min).  Fig.  5  further  shows  several  distinct  features 
which  can  be  easily  identified.  Firstly,  the  stress  in  the  electrolyte 
and  the  cathode  are  always  compressive  while  there  are  both 
compressive  and  tensile  stresses  in  the  anode.  Secondly,  both  the 
maximum  compressive  stress  across  the  three  layers  and  the 
maximum  tensile  stress  in  the  anode  layer  increase  as  the  heat-up 
time  decreases.  The  reason  is  that  the  temperature  gradient 
through  the  fuel  cell  plate  increases  significantly  when  the  heat-up 
time  decreases,  leading  to  greater  thermal  stress  to  the  cell. 

It  is  a  well-known  fact  that  the  ceramic  materials  are  susceptible 
to  fracture  under  tensile  stress  and  the  compressive  fracture 
strength  of  a  ceramic  material  is  normally  much  greater  than  its 
tensile  fracture  strength.  The  compressive  stress  in  the  cathode  is 
kept  in  a  low  level  so  the  cathode  is  usually  unlikely  to  crack.  The 
compressive  stress  within  the  electrolyte  layer  is  much  higher, 
which  can  even  reach  800  MPa  although  the  stress  could  be  kept 
below  the  critical  compressive  fracture  strength  of  8YSZ  which  has 
been  reported  to  be  greater  than  1  GPa  [30].  Therefore,  the  anode 
layer  which  experiences  the  tensile  stress  is  most  susceptible  to 
fracture  in  this  situation. 

The  failure  probability  of  the  fuel  cell  can  be  calculated  by 
Weibull  statistics  [31,32  as  following: 


Pel-  exp 


(31) 


where  ao  and  m  are  the  characteristic  strength  and  Weibull 
modulus,  respectively.  The  fracture  strength  properties  of  the  ma¬ 
terials  at  operation  temperature  are  summarized  in  Table  8. 


Table  9 

Maximum  tensile  stress  and  failure  probability  of  different  heat-up  rates. 


DFFC  case 

SOFC  case 

Heat-up  time  (s) 

0.1 

1 

50 

1000 

(Tmax  in  anode  (MPa) 

83 

18 

10.8 

8.7 

Failure  probability 

0.76 

3.2E-05 

8.9E-07 

2.0E-07 

The  maximum  tensile  stress  and  failure  probability  are  shown  in 
Table  9.  The  rapid  heat-up  in  DFFC  operation  may  lead  to  large 
stress  to  the  cell  due  to  the  temperature  difference  through  the  cell 
thickness.  The  corresponding  failure  probability  shows  four  times 
of  magnitude  increase,  from  2.0E-07  to  8.9E-07,  when  the  operation 
condition  transfers  from  “SOFC  case”  to  “DFFC”  case,  i.e.  the  heat-up 
time  decreases  from  1000  s  to  50  s.  And  it  shows  two  orders  of 
magnitude  increase  when  the  heat-up  time  is  1  s.  What’s  more, 
there  exists  another  four  orders  of  magnitude  increase  when  the 
cell  experiences  a  faster  heat-up  rate,  which  is  0.1  s  from  298  K  to 
1073  K.  As  a  conclusion,  the  rapid  heat-up  rate  in  DFFC  case  lead  to 
larger  thermal  stress  to  the  cell,  thus  increasing  the  failure  proba¬ 
bility.  However,  by  controlling  the  heat-up  rate  at  a  relatively  slow 
level,  the  increase  of  stress  and  failure  probability  can  fall  into  an 
acceptable  range. 


3.3.  Stress  comparison  for  different  fuel  cell  configurations 

From  the  previous  discussions,  it  can  be  seen  that  the  cell  in 
DFFC  operation  is  more  susceptible  to  fracture  compared  to  the  cell 
for  being  used  in  common  SOFC  operating  condition  due  to  its  rapid 
heat-up  rate.  It  is  therefore  valuable  to  choose  an  appropriate  SOFC 
configuration  which  is  more  suitable  for  DFFC  from  the  point  of 
thermal  shock  resistance.  In  the  present  study,  two  different  SOFC 
configurations,  the  ASSOFC  and  the  ESSOFC,  are  compared.  In  the 
model,  they  get  through  the  same  heat-up  process  within  1  s  from 
298  K  to  1073  I<  and  the  temporary  temperature  and  stress  are 
calculated  when  the  anode  reaches  the  operation  temperature.  The 
transient  temperature  distribution  and  stress  distribution  are 
shown  in  Figs.  6  and  7. 

When  the  anode  temperature  reaches  the  usual  SOFC  operation 
temperature  1073  K,  the  transient  temperature  difference  between 
the  cathode  and  the  anode  is  about  20  K  in  ESSOFC,  which  is  much 
smaller  than  that  in  ASSOFC  (more  than  100  I<).  The  reason  is  that 
the  cell  thickness  of  ESSOFC  configuration  is  much  smaller  than 
that  of  ASSOFC  and  this  contributes  to  the  rapid  heat  conduction  in 
the  thickness  direction  and  also  helps  to  decrease  the  temperature 
difference  between  the  cell  layers.  Despite  of  the  smaller  temper¬ 
ature  difference,  the  maximum  tensile  stress  in  ESSOFC  is  greater 
than  that  in  ASSOFC.  And  the  cathode  layer  of  ESSOFC  undergoes  a 
tensile  stress  which  is  different  from  the  situation  in  ASSOFC,  which 
means  that  the  cathode  of  ESSOFC  is  more  susceptible  to  fracture 
compared  to  ASSOFC.  The  failure  probability  calculation  results  of 
both  cases  are  shown  in  Table  10.  The  failure  probability  of  ESSOFC 
magnifies  two  orders  compared  to  ASSOFC  when  undergoing  a 
sudden  heat-up  process,  which  means  that  ASSOFC  has  better 
thermal  shock  resistance  than  ESSOFC  thus  is  more  suitable  in  DFFC 
operation. 


3.4.  Failure  probability  analysis  for  DFFC  operation  period 

Except  for  the  differences  between  operating  conditions  of  DFFC 
and  common  SOFC  due  to  the  rather  rapid  heat-up  rate,  another 
significant  difference  is  the  operating  environment  including  both 
the  temperature  and  the  composition.  For  example,  the  operating 
temperature  of  the  common  SOFC  is  well  controlled  by  the  electric 
furnace  while  the  temperature  of  DFFC  is  mainly  determined  by  the 
flame.  In  fact,  the  flame  temperature  and  composition  may  vary 
significantly  in  a  plane,  such  as  the  case  of  the  flame  provided  by  an 
ethanol  lamp  [2  or  a  Bunsen  burner  [1  ].  In  this  section  the  effect  of 
flame  non-uniformity  on  both  the  performance  and  the  stress  filed 
of  the  fuel  cell  is  studied.  The  temperature  distribution  is  supposed 
to  follow  the  same  trend  as  shown  in  Fig.  8.  The  degree  of  the  non¬ 
uniformity  is  characterized  by  the  various  temperature  difference 
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Fig.  6.  Transient  temperature  distribution  (a)  ASSOFC  (b)  ESSOFC. 


Fig.  7.  Transient  stress  distribution  (a)  ASSOFC  (b)  ESSOFC 


Table  10 

Maximum  tensile  stress  and  failure  probability  for  ASSOFC  and  ESSOFC. 


ASSOFC 

ESSOFC 

Anode 

Anode 

Cathode 

ffmax  (MPa) 

18 

35 

20 

Failure  probability 

3.2E-05 

3.3E-03 

7.5E-03 

Fig.  8.  Temperature  distribution  in  x  direction. 


AT.  The  anode  gas  composition  is  set  as  the  equilibrium  flame 
composition  at  0  =  1.4,  as  shown  in  Table  11. 

Table  12  shows  the  maximum  tensile  stress  and  the  corre¬ 
sponding  failure  probability  of  various  AT  and  Fig.  9  shows  the  ionic 
current  density  distribution  of  various  AT  at  a  same  voltage  0.7  V. 
The  flame  non-uniformity  induces  significant  thermal  stress  to  the 
fuel  cell  and  increases  the  failure  probability  greatly.  It  should  be 
noted  that  when  the  AT  reaches  750  K,  i.e.  the  flame  does  not  cover 
the  whole  fuel  cell  area,  the  failure  probability  reaches  1.  It  should 
be  noted  that  the  fuel  cell  performance  is  also  greatly  influenced  by 
the  flame  temperature  distribution,  which  is  shown  in  Fig.  9.  It  can 
be  seen  that  the  ionic  current  density  increases  rapidly  near  the 
electrolyte,  which  means  that  the  electrochemical  reaction  rate 
grows  quickly  near  the  electrolyte.  Fig.  9  also  indicates  that  the 
distribution  of  the  current  density  is  non-uniform,  which  is 
consistent  with  the  uneven  flame  temperature  distribution.  The 
maximum  ionic  current  density  exists  where  the  flame  tempera¬ 
ture  reaches  the  highest  and  decreases  rapidly  in  x  direction  as  the 
temperature  decreases.  And  the  degree  of  the  ionic  current  density 
non-uniformity  is  dependent  on  the  flame  temperature  uneven 
degree. 


Table  11 

Equilibrium  flame  composition  at  </>  =  1.4 


h2  (%) 

h2o  (%) 

CO  (%) 

C02 (%) 

n2  (%) 

6.0 

17.5 

7.5 

4.8 

64.2 

Table  12 

Maximum  tensile  stress  and  failure  probability  of  various  AT 


AT(K) 

100 

200 

750 

O'max  (MPa) 

38 

56 

196 

Failure  probability 

5.9E-03 

8.6E-02 

1 
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Fig.  9.  Current  density  distribution  of  various  AT  (a)  AT  =  100  K  (b)  AT  =  200  K  (c) 
AT  =  750  K. 


3.5.  Failure  probability  analysis  for  DFFC  start-up  period 


Fig.  11.  Transient  stress  field  in  the  DFFC  start-up  period. 


that  the  transient  temperature  shows  great  non-uniformity  in  the 
fuel  cell,  in  both  the  x  and  y  directions,  bringing  huge  stress  to  the 
cell.  The  failure  probability  reaches  8E-02,  which  is  much  bigger 
than  the  situation  that  only  one  factor  is  taken  into  account  (the 
rapid  heat-up  rate  or  the  temperature  non-uniformity).  From  the 
results  above,  we  can  come  to  a  conclusion  that  the  start-up  period 
is  the  time  when  the  cell  is  most  likely  to  fracture,  and  the  flame 
temperature  uniformity  is  rather  vital  in  DFFC  operation,  especially 
in  the  start-up  period. 


4.  Conclusion 

A  detailed  two-dimensional  model  of  direct  flame  fuel  cell  was 
developed  by  considering  the  coupling  effects  of  heterogeneous 
chemical  and  electrochemical  reactions,  electrode  microstructure, 
transport  processes  of  mass  and  charge,  as  well  as  the  thermal 
mechanical  stress.  The  stress  distribution  was  simulated  at 
different  heat-up  rates  which  represents  the  typical  DFFC  and  the 
common  SOFC  operation.  Both  the  ASSOFC  and  the  ESSOFC  are 
simulated  under  the  same  heat-up  condition  for  DFFC.  The  failure 
probability  for  both  the  DFFC  operation  period  and  the  start-up 
period  was  analyzed.  The  model  is  demonstrated  to  be  a  useful 
tool  for  understanding  the  mechanical  stress  distribution  within 
DFFC  cell  and  for  the  cell  structure  design  and  optimization.  The 
results  indicate  that: 


In  the  start-up  period  of  DFFC,  the  fuel  cell  will  not  only  get 
through  a  rapid  heat-up  process,  the  flame  temperature  can  also  be 
uneven.  In  this  section,  the  heat-up  rate  is  set  to  be  1  s  from  room 
temperature  to  operating  temperature  and  the  temperature  dif¬ 
ference  AT  defined  in  Section  3.4  is  set  100  K.  The  transient  tem¬ 
perature  and  stress  field  are  shown  in  Figs.  10  and  11.  It  can  be  seen 


(1)  The  rapid  heat-up  process  which  exists  in  DFFC  operation 
brought  in  a  great  temperature  gradient  in  the  cell  thickness 
direction,  increasing  the  failure  probability  by  2-6  orders. 
But  by  controlling  the  heat-up  rate  at  a  relatively  slow  level 
(50  s  from  298  K  to  1073  I<),  the  increase  of  stress  and  failure 
probability  can  fall  into  an  acceptable  range. 

(2)  The  failure  probability  caused  by  thermal  shock  of  ESSOFC  is 
about  two  orders  higher  than  that  of  ASSOFC  for  DFFC 
operation.  It  suggests  that  ASSOFC  would  be  better  for  being 
utilized  in  DFFC  case  in  the  aspect  of  thermal  shock 
resistance. 

(3)  The  non-uniform  flame  temperature  will  not  only  induce 
more  significant  thermal  stress,  but  also  lead  to  serious  non- 
uniform  cell  performance  which  may  greatly  increase  the  cell 
failure  probability. 

(4)  In  the  start-up  period,  the  non-uniform  temperature  distri¬ 
bution  of  the  cell  due  to  the  combined  effects  of  the  quickly 
rising  of  the  temperature  and  the  non-uniform  flame  tem¬ 
perature  exists  both  in  the  cell  thickness  direction  and  the 
cell  plane  direction,  further  increasing  the  failure  probability 
to  around  8E-02.  This  further  indicates  a  flame  which  could 
provide  a  uniform  temperature  field  is  of  great  importance, 
especially  during  the  start-up  period. 
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